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Abstract In the following we consider a 2-dimensional 
system of ODE's containing quasiperiodic terms. The 
system is proposed as an extension of Mathieu-type 
equations to higher dimensions, with emphasis on how 
resonance between the internal frequencies leads to a 
loss of stability. The 2-d system has two 'natural' fre- 
quencies when the time dependent terms are switched 
off, and it is internally driven by quasiperiodic terms in 
the same frequencies. Stability charts in the parameter 
space are generated first using numerical simulations 
and Floquet theory. While some instability regions are 
easy to anticipate, there are some surprises: within in- 
stability zones small islands of stability develop, and 
unusual 'arcs' of instability arise also. The transition 
curves are analyzed using the method of harmonic bal- 
ance, and we find we can use this method to easily 
predict the 'resonance curves' from which bands of in- 
stability emanate. In addition, the method of multiple 
scales is used to examine the islands of stability near 
the 1:1 resonance. 

Keywords 2-d Mathieu • quasiperiodic • Floquet • 
harmonic balance • multiple scales 



1 Introduction 

The study of systems of ode's with periodic coefficients 
arises naturally in many branches of applied mathemat- 
ics. Aside from being interesting in their own right, in 
applications such systems arise when considering the or- 
bital stability of a periodic solution to some dynamical 
system (see, for example, Betounes [1]). The archetypal 
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equation for this class is the Mathieu equation, 
X + {6 + ecost)x = 



(1) 



where a dot denotes differentiation w.r.t. time. This 
equation, though linear, is sufficiently complicated by 
its explicit time dependence that closed form solutions 
do not exist. Instead, the common method of attack is 
to consider the stability of the origin using Floquet the- 
ory (or other methods) . The {6, e) parameter space is 
clearly divided into regions of stability/instability (see 
Jordan and Smith [5] for a detailed discussion) . The ap- 
pearance of these instability zones is intuitively under- 
stood in the following way: consider the time-dependent 
term in (HI) as a perturbation. When e = 0, x has pe- 
riodic solutions with frequency %A5 (the 'natural' fre- 
quency of the system). As we switch on the pertur- 
bation, there will be a resonance between the 'driving' 
frequency, 1, and the natural frequency vS, which leads 
to instability. 

There are many extensions to this equation, of which 
we describe only a few. Firstly the equation can be 
extended to many dimensions, by replacing x with a 
vector and S and e with matrices. This problem was 
analyzed in Hansen [3] using Floquet theory, however 
(in contrast with the present work) the time dependent 
terms only depended on one frequency. Secondly the 
equation may be extended to contain nonlinear terms; 
for example El-Dib 2,, consider an extension up to cu- 
bic order with each term containing subharmonics in 
the periodic term, and Younesian et al [13] consider 
also cubic order with two small parameters and terms 
in sin and cos. 

Thirdly, and most importantly for the present work, 
we may develop the time dependent term to contain 
two periodic terms in different frequencies, and as the 
two frequencies need not be commensurable this system 
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is referred to as the 'quasiperiodic Mathieu equation'. 
This type of system was developed and studied in great 
detail by Rand and coworkers: in Rand and Zounes [M] 
and Rand et al [8] the following quasiperiodic Mathieu 
equation is examined 

x + [5 + e{cos{t) + cos{ujt))]x = (2) 

and transition curves in the {S, uj) parameter space are 
developed. Specific resonances are examined in Rand 
and Morrison [7] and Rand et al |Bj. In Sah et al [S] a 
form of the quasiperiodic Mathieu equation is derived 
and examined in the context of the stability of motion 
of a particle constrained by a set of linear springs, and 
finally in Rand and Zounes |15j a nonlinear version of 
the quasiperiodic Mathieu equation is studied. In these 
papers a combination of numerical and analytical tech- 
niques are used to examine how resonances between the 
internal driving frequencies lead to a loss of stability in 
solutions, and it is very much in the spirit of these pa- 
pers that we present the following work. 

We seek to extend consideration of quasiperiodic 
Mathieu- type systems to higher dimensions. The mo- 
tivation for this is clear: as mentioned previously, lin- 
ear systems with (quasi)periodic coefficients often arise 
when a nonlinear dynamical system is linearized around 
a (quasi)periodic solution, the so-called 'variational equa- 
tions' [1]. As most dynamical systems of interest will 
have more than one degree of freedom (for example the 
motion of a particle in space under the influence of a col- 
lection of forces), the associated variational equations 
will also have more than one degree of freedom (unless 
only certain modes are being considered as in [9]). We 
seek a 2-dimensional system containing more than one 
frequency and a small parameter e to enable a pertur- 
bation analysis. However, it is beneficial to keep the 
number of free parameters low in order to make the be- 
havior clear. With this in mind, we shall consider the 
following system: 

X + [a'^ + e cos{(3t)]x + e cos{at)y — 
y + e cos{f3t)x + [(3'^ + e cos(at)]y 
which in general is of the form 



(3) 



[Po + ePiit)]X^O 



(4) 



^0 



(5) 



with X — (x, y) and 

a^ \ p (i\ — f cos{/3t) cos(ai)" 
/3V ' ^1^*^- [cosipt) cos{at) ^ 
Here a dot denotes differentiation w.r.t. time, e is small 
and a, /3 G R. 

We can predict to some degree how solutions of this 
system will behave: when e — 0, the solutions are peri- 
odic in the natural frequencies a and (3. As the pertur- 
bation, ePi{t)X, is switched on we would expect reso- 
nance between the natural and driving frequencies, as 



in the Mathieu equation (IT|) . As the natural and driving 
frequencies are the same, this means we expect insta- 
bility to arise where a and /3 are commensurable. And 
this is to a degree what we do see, with two unexpected 
features. Firstly, curved 'arcs' of instability can be seen 
which clearly do not arise from integer ratios of the pa- 
rameters. Secondly, within the bands of instability small 
islands of stability appear. These two features will be 
addressed in more detail in what follows. 

While this system may seem a little artificial, we will 
show how it may arise as an approximate variational 
equation. Consider the following system of nonlinear 
ODE's: 



-a u 



-P^v- 



uv. 



(6) 



This is an example of a system which is not dissipa- 
tive but nonetheless is not Hamiltonian; as such it falls 
into the more general class of Louiville systems (see be- 
low). The author's direct experience with systems of 
this type is in the problem of the solar sail in the cir- 
cular restricted 3-body problem (see for example the 
author and Mclnnes [TOj), although other applications 
exist. 

The origin of ([5]) is a fixed point, and letting u = 
u + x, V — V + y we have the linear system 

d^X 



dt^ 






X 



(7) 



u^u,v^v 



where X — {x,y). Letting (u,w) — (0,0) we find two 
linear oscillators which we write as ecos(ai), ecos{/3t) 
(letting the amplitudes be equal and setting the phases 
to zero), and then letting u — ecos{at), v — ecos{/3t) 
we get system ([3]) above. Note that the autonomous 
system ([6]) is not Hamiltonian and therefore (O is not 
symmetric; neither, for the same reason, is Pi{t) in 
([3]). Whether or not the system under consideration is 
Hamiltonian only becomes an issue in dimension 2 and 
above; in one dimension we can always find a (time 
dependent) Hamiltonian, a fact taken advantage of in 
Ye§ilta§ and §im§ek [T^ . 

Note also that {u,v) = (ecos(ai),ecos(/3i)) is not a 
solution to ([ni, rather a solution to the linearized sys- 
tem (O. As such ([3]) is not a true variational equation, 
rather an approximation to one. We make this approxi- 
mation, rather than seek more accurate solutions to ([B]). 
so as to enable an analytical rather than purely numer- 
ical treatment. The Mathieu equation itself can be seen 
in this light, as an approximate variational equation to 
the system 



The fact that ^ is only approximately a variational 
equation will be of relevance in what follows. 
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The rest of the paper is laid out as follows: In Sec- 
tion 2 we numerically integrate system ^ to examine 
the stability properties of various values of a and /?. 
As computers can only handle finite precision we must 
choose rational values of a and /3 and therefore ^ is 
periodic, rather than quasiperiodic (a more appropriate 
term might be multiply-periodic). As the system is pe- 
riodic, albeit of possibly very long period, we may use 
Floquet theory to enhance the numerical integration. 
We are able to generate stability charts in the (a, /3) pa- 
rameter space for various values of e and we discuss the 
main features of same. The fact that we have lost the 
true quasiperiodic nature of the system in this section 
is made up for by the following: every irrational number 
has rational numbers arbitrarily close to it which thus 
give a very good approximation; the numerically gener- 
ated stability charts allow us to isolate the interesting 
features of the system which informs later analysis; and 
the remaining analysis in the paper holds for arbitrary 
a and (3 and thus recaptures the quasiperiodic nature 
of the system. 



In Section 3 we use the method of harmonic balance 
to find the transition curves in parameter space. Writ- 
ing solutions in terms of truncated Fourier series we 
approximate the transition curves by the vanishing of 
the Hill's determinant. As the truncation size becomes 
large the matrices of vanishing determinant also become 
large making this method cumbersome. Also, there are 
convergence issues with the infinite Hill's determinant 
appropriate to ^ which we will elaborate on. How- 
ever, we show that we can use the method of harmonic 
balance to quickly find approximations to transition 
curves by understanding the following mechanism of 
instability: when e = there are certain 'resonance 
curves' in parameter space (for the Mathieu equation 
([T]) these are simply the points 5 = n^/4 for n G Z). 
As £ grows, these curves widen into bands of instabil- 
ity. The method of harmonic balance allows us to find 
these resonance curves quickly, and we demonstrate the 
usefulness of this technique on a second system closely 
related to ([3]). 



In Section 4 we use the method of multiple scales 
to examine in detail the fine structure near the a ~ /3 
resonance. From Section 1 we can see small islands of 
stability arising in the middle of a large band of in- 
stability, and we examine this feature more closely. We 
separate out the 'slow flow' using three time variables, 
and find we can analytically predict the appearance of 
these pockets of stability. Finally we give some discus- 
sion and suggestions for future work. 



2 Numerical integration 

We wish to fix the values of the parameters {a,f3,e) 
and numerically integrate the system ^ to examine 
the stability of the origin. In deciding what range of 
values to choose for the parameters, we note that Q 
has the following scale invariance: consider the system 

X + [{ma) + m e cos{m/3t)]x + m e cos(mat)y — 

y + rn?e cos{mf5t)x + [(to/3)^ + rn?e cos{mat)]y — (8) 

By defining a new time coordinate t' = mt we recover 
the original system ^. Thus the following parameter 
values are equivalent as regards stability: 

(a, /3,e), {ma, 171(3,171^8). 

This means we need only consider a and /3 on the unit 
square (for example) , and the stability of other values 
of a and /3 are found by varying the value of e. For 
example, the stability of system ([3]) when a — 3, /3 = 
2,£ = 0.1 is the same as when a = 3/4, /3 = 1/2, e = 
0.0125. 

As mentioned in the introduction, to integrate nu- 
merically we must use finite precision values of a and 
13 and thus the system is (multiply-)periodic. As such, 
we may take advantage of Floquet theory. 



Writing ([3|) as a first order system, the fundamental 
matrix ${t) solves the following initial value problem 
(see, for example, Yakubovich and Starzhinskii [11] for 
a detailed treatment of Floquet theory) 
d<P{t) 



= A{t)<P{t), ^{0)^1 



(9) 



A{t) = 



(10) 



dt 
where 

/ 10\ 

1 

— [a^ + £Cos(/3t)] — £cos(Q;t) 
V -£cos(/3i) -[(3^ + e cos{at)] / 

and A{t + T) = A{t) for some fundamental period T. 
The stability of the solutions to this system are deter- 
mined by the monodromy matrix|j, given by 'P{T), since 
it can be shown that 



x{T) 

y{T) 



^ '[yio) 



More precisely, stability is determined by the eigenval- 
ues of 'P{T), the so-called (characteristic) multipliers 
A;. If I Ail < 1 the solutions remain bounded, implying 
stability, and if any |Ai| > 1 the solution is unbounded. 
From the structure of (jlOp we can make some com- 
ments about the multipliers: 

^ fn some of the literature, A{t) is known as the monodromy 
matrix and <P{T) as the matrizant, propagator or character- 
istic matrix. 
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(a) (6) (c) id) 

Fig. 1 Possible positions of multipliers in the complex plane with respect to the unit circle for a t-invariant system of order 4. 



1. Louiville's theorem [5] gives 

det[^{T)] = A1A2A3A4 = exp j / Tr{A{t)}dt] = 1. 

This however does not mean that the eigenvalues 
must appear in reciprocal pairs, merely that all four 
multiplied together must be 1. 

2. That they must appear in reciprocal pairs comes 
from the t-invariant nature of the system TTD . These 
are systems with the property 



A{-t)G + GA{t) = 



(11) 



for some non-singular constant matrix G. Since A{t) 
is even, i.e. A{—t) — A{t), this holds with 



G = 



I 
0-/ 



(12) 



the entries understood to be 2 x 2 matrices. In this 
case, we can show 



<P{T)-^ = G.<P{T).G- 



(13) 



in other words, the monodromy matrix and its in- 
verse are similar and thus have the same eigenvalues; 
thus if Xi is a multiplier then so is l/A^. 
As ^{T) is real, complex eigenvalues must appear in 
conjugate pairs, and since they must also be recip- 
rocals this means complex eigenvalues must be on 
the unit circle (with one rare exception, see point 
4 below). Therefore, complex multipliers represent 
stability, and real multipliers instability. 
If the system ^ were a variational equation, that is 
a (nonlinear autonomous) system linearized around 
a periodic solution, we would expect a unit eigen- 
value (see chap. 7 of Betounes [T]) and therefore 
(from the previous point) 2 unit eigenvalues. This 
would make life easier, as stability would be ensured 
if 

Tr{^{T)} < 4 (14) 

and thus we would not need to calculate the multi- 
pliers themselves. Unfortunately we cannot say this 



is the case (see the Introduction) , and so there is no 
guarantee of unit eigenvalues. As the system can be 
viewed as an approximation to a variational equa- 
tion then we would expect (for small e) to have two 
multipliers close to, but not equal to, one. 
4. There are a number of possibilities for the positions 
of the multipliers in the complex plane, shown in 
Figured] The most common is for the set 

1 + e, , a + bi, a — bi (15) 

1 + e ^ ' 

with e e R+ and a^ -t- 6^ ^ 1 (since A = 1/A). We 
may describe this as 'saddle-centre' (Figure [Ub)). 
Stability is given by 'centre-centre' configurations 
(Figure [He)), however we must also accept the pos- 
sibility that the multipliers may not be on the unit 
circle (Figure [ijd)). This is possible when the eigen- 
values of the centre-centre configuration come to- 
gether and then bifurcate off the unit circle (a Krein 
collision). Note, the intermediate stage here means 
the eigenvalues are equal and this is a degenerate 
case; this will not be pursued further. 

Based on these considerations, we will calculate the 
multipliers of the monodromy matrix and find their 
norms; if the norm is less than some cut-off value we 
say the system is stable, otherwise unstable. 

We focus on (a, /3) = (0, 1) x (0, 1) (as we may do due 
to the scale invariance of the system, mentioned previ- 
ously). Splitting this interval into n equal segments, we 
let a = i/n and (3 = j /n for i,j = 1, . . . ,n. The funda- 
mental period of system (fTUjl is given by 



T: 



27rn 



GCD[z,j] 



where GCD is the greatest common divisor. We then 
integrate the system (jH]) for T units of time, evaluate 
<?(r), and calculate its eigenvalues. 

The most important issue with regards accuracy is 
the long integration time. If i and j do not have a com- 
mon divisor, than T = 27rn, which is large for n large. 
We may halve this integration time by taking advantage 
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Fig. 2 Stability diagram for system JSj with £ = 0.1 and the a, /3 unit square divided into an 800 x 800 grid. White represents 
points where the norm of the maximum multiplier is greater than 1.025 (and therefore the system is unstable), and black 
where the norm is less than 1.025 (the system being stable). The dashed circle shows the region of parameter space where the 
solutions in Figure [3] are taken from. 





Fig. 3 Long-time solutions for very close together points in a, 13 parameter space. On the left is the solution for a = 
435/800, /3 = 425/800; this point is in the stability island indicated in Figure[2] On the right is the solution for a = 436/800, /3 = 
425/800; this solution is outside the stability island. The integration times here are IBOOtt, and in both columns we present the 
X- and j/-solution versus time and the x-y solution. The initial conditions for both are x{0) = l,y(0) = l,i(0) = 0,y{Q) = 0. 
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of the t-invariant nature of the system (see above). If 
A{t) is such that dTTl) holds, then we may show that 

<?(r) = G.^{T/2y^.G.^T/2) 

since G = G^^. Hence we need only integrate over one 
half period to calculate 'P{T). This improves the accu- 
racy greatly, with one caveat: we must invert <P[T /2). 
However, since Louiville's theorem holds for all t, this 
means the determinant of (p{T /2) should be one, and 
thus inverting may be done accurately. 

In fact, we may use the condition det[<P{t)\ = 1 to 
monitor the accuracy of the integration, however this 
proves very costly with regards integration time. In- 
deed, for some values of a, /? and e we would expect 
the solution to be grossly unstable. Monitoring the ac- 
curacy of such solutions is inefficient as our concern is 
the stability-instability divide, rather than how large 
these very large multipliers may be. As such, we are 
satisfied that if det[(l>{t)] grows very large this is due to 
the huge growth in very unstable combinations of a and 
/3; we will not waste effort with the accuracy of these 
solutions. 

In Figure [5] we give the stability diagram for system 
([31) with £ = 0.1. Here, we set n — 800 and thus divide 
the a, P unit square parameter space into 800^ data 
points. For each we calculate the maximum multiplier 
norm, with white regions being where the max norm is 
larger than 1.025, black regions where the max norm is 
less than 1.025. We may make some observations: 

— The most striking features are the large, linear bands 
of instability. These are aligned along lines of con- 
stant slope in the parameter space, in other words 
where a//3 =constant. These are clearly the domi- 
nant, low order resonances between a and /3, and the 
lower the order of resonance the wider the band of 
instability. There are also higher order resonances, 
and these are of narrower width. This behavior is 
very much what we would expect, based on the clas- 
sical Mathieu equation and its variations. 

— There are two interesting and unexpected features. 
The first is the large 'arc' of instability seen in the 
corners. This curved band of instability is not pre- 
dicted by simple resonance between a and /3. 

— Second and perhaps more interesting are the small 
islands of stability which appear in the large a ~ /? 
resonance band. We would fully expect the bands of 
instability to grow with e, and this we do observe. 
However, for small pockets of stability to develop 
when surrounding them is very large regions of in- 
stability is most unexpected. These islands are pre- 
carious: we show in Figure [3] the solution for two 
very nearby points in parameter space, one inside a 
small island of stability and one outside. 



— The corresponding plot for other values of e can be 
inferred from Figure [2] in the following way. Larger 
values of e generate stability charts which 'zoom 
in' towards the origin of Figure [2l smaller values 
'zoom out'; this is due to the scale invariance of the 
system as described at the beginning of this section. 
For example, the stability chart for a, /? on the unit 
square and e ~ 0.4 is the same as the lower left 
quadrant of the same chart for e = 0.1, since the 
parameter values (a, /3, e) = (1/2, 1/2, 0.1) have the 
same stability properties as (1,1,0.4). This is why 
there is a cluster of unstable points near the origin: 
moving towards the origin is equivalent to increasing 
£; the resonance bands widen and overlap. As we 
decrease e, regions beyond the unit square in Figure 
[2] are drawn in; since the resonance bands narrow as 
we move away from the origin this is equivalent to 
saying they narrow as e is decreased, as expected, 
and the chart is increasingly made up primarily of 
stable points. The exception is the a ~ /S resonance 
band; it does not appear to taper as we move along 
it away from the origin. Some insight into this is 
gained in Section 4. 

While the plots generated using numerical integra- 
tion are informative, they give us no feel for why certain 
bands of instability develop or which are dominant. Nor 
do they give us an analytic expression for the values at 
which the system bifurcates from stable to unstable. 
In the next two sections, we use analytical techniques 
which do just this, and what's more recapture the true 
quasiperiodic nature of the original system. 



3 Transition curves using the method of 
harmonic balance 

In Floquet theory, t-invariant systems must have multi- 
pliers on the unit circle or real line (with the rare excep- 
tion of Figure [T](d)). As the system passes from stable 
to unstable, the multipliers must pass through the point 
(1,0) or (—1,0) in the complex plane. These multipli- 
ers represent solutions which have the same period, or 
twice the period, of the original system. This means so- 
lutions transiting from stable to unstable have a known 
frequency and thus can be written as a Fourier series. 
For this series to be a solution to the system puts a 
constraint on the parameters involved, and this in turn 
tells us the transition curves in parameter space. 

In quasiperiodic systems, such as the one under con- 
sideration in this paper, which we give again as 



X +[a'^ + e cos{f5t)]x + e cos{at)y = 
y + e cos{l3t)x -f [/3^ -I- e cos{at)]y = 0, 



(16) 
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Floquet theory does not hold. However, based on the 
discussion above we can use the following ansatz: solu- 
tions on transition curves from regions of stability to 
instability in parameter space have the same quasiperi- 
odic nature as the system itself. This was justified and 
used by Rand et al in [H1I5] , and while lacking in rigor- 
ous theoretical proof we find this assertion works well 
in predicting the transition curves in parameter space. 
Since any quasiperiodic function can be written as 
an infinite Fourier series in its frequencies (see, for ex- 
ample, Goldstein et al [3 ), we write the solutions to 
system ^TE\\ in the following form, 



E 



-oo 7n^ — oo 



exp 



— [na + mp) 



(17) 



where the 2 appears in the exponent to account for 
solutions with twice the period of the system (corre- 
sponding to the —1 multiplier). We substitute (fT7|) into 



TB|) . use the identity cos{9t) = h{e 



i9t 



i9t 



) , and shift 



indices to cancel out the exponential terms. What re- 
mains are the following two recurrence relations: 



{na + mpy 



Ar, 



+ 9 (^ri,m-2 + An^m+2 + Bn-2,m + -B„+2,m) — 0, 



(18) 



2 [na + mPY 



Br, 



+ 9 {^n,m-2 + ^n,m+2 + Bn-2,m + -Bn+2,m) — 0, 

(19) 

with n,m = — cx), . . . , oo. This is a set of linear, homo- 
geneous equations in Anm , Bnm of the form 

Cx = 0, 

and for non-trivial solutions to exist we require det{C) — 
0; these are the well know infinite Hill's determinants. 
In practice, we truncate at some order N so as n^m — 
—N,...,N. Truncating at larger values of N will pro- 
duce more accurate determinants, provided the deter- 
minant converges. However, close inspection of the re- 
currence relations (|18I19|) reveals that this is an issue 
for the system in question. 

Writing the coefficient of Anm in (dH as jnm and 
dividing across gives the following: 



Ar 



27r. 



{An.r> 



= 



(20) 



with a similar expression for (J19p . In most systems, as 
TV — > oo the 'off-diagonal' terms here grow progressively 




Fig. 4 Some transition curves on the unit square for e = 0.1 
using the method of harmonic balance. The curves are where 
the determinant of the coefficient matrix C vanishes. 



smaller, ensuring convergence. However, we see that if 
we let, for example, n — 2,m = then 

720 = 0. 

Thus the infinite determinant does not converge. This 
is a problem particular to systems with a driving fre- 
quency that is equal to a natural frequency. 

As a way around this problem we notice that a so- 
lution with period 2T also has period 4T, and indeed 
6r, ST, . . . and so on. Thus we may write 



E 



E 



n— — oo m— — oo 



exp 



it 
2M 



(na + mj3) 



(21) 



where M = 1,2,3,.. .. In this case we see 



In 



(na + 771/3)^ 
(2M)2 



and the singular term does not appear in the Hill de- 
terminant until N — 2M. Thus we may work up to 
order N = 2M — 1 and calculate the determinant. It 
will be a function of the parameters {a,/3,e) and by 
setting it equal to zero we may therefore find analytic 
approximations to the transition curves in parameter 
space. 

For example, with M = 3 we generate the matrix 
appropriate to 5th order, find its determinant, and set 
it equal to zero. Factorizing the determinant we see it 
contains terms like the following: 

(4a2 ~(3^ + e){Aa^ - /^^ _ g) ^ o, 
(9a* - 82a2;32 ^ 9^4 _ ^^^2-^ ^ 

(9a* - 82a2/32 + 9/3* - IGe^) = 0. 
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The first two clearly indicate a band of instability with 
width 0{s) around the 2:1 resonance curve, which grows 
to either side of the resonance as e grows. The second 
two terms indicate a band of width O(e^) around the 3:1 
resonance curve, which grows away from the resonance 
as e grows. This implies that stability for resonant val- 
ues of a and /? is not as straightforward as we would 
initially expect: for a/P — a/b with a, b integers the 
band of instability may move away from the resonance 
curve as e grows and thus commensurable values oia, (3 
will in fact be stable. 

We give the zeroes of some of the main factors of 
the determinant in Figure HI We can see in the figure 
some of the main features we expect from our numerical 
simulations, however some features are notable absent; 
foremost among these is the 1:1 resonance. Going to 
higher orders will capture more of the transition curves, 
however the dimension of the problem quickly becomes 
very large. If we truncate the series at iV = 2M — 1 
we must find the determinant of a 2{2N + 1)^ square 
matrix, which for the example described above (M = 
3) is 242 X 242. Calculating the determinant can be 
facilitated in the following way: these large matrices 
are very sparse as is usual for recurrence relations. We 
may separate out 2M submatrices, each describing the 
T,2T,..., 2MT solutions, as in 

C = C1C2 . . . C2M, 

and set the determinants of each to zero. However, the 
improvements in accuracy from increasing M do not 
justify the increasingly cumbersome and lengthy ma- 
trices. Thus for this system it is worthwhile to study 
the resonances of greatest import individually; this we 
do in the next section. 



Aside from these complications, we note that the 
method of harmonic balance can be used very quickly 
and directly to predict, in broad terms, the locations 
of the bands of instability. As mentioned previously, we 
may understand the mechanism of instability so: reso- 
nance curves in the (a, /3) parameter space grow with e 
into bands of instability due to resonance between the 
natural and driving frequencies. Predicting the location 
of said resonance curves will give approximately the re- 
gions of parameter space which will become unstable. 
As these resonance curves are the bands of instability 
in the limit of e — >■ 0, we may simply let e = in the 
recurrence relations given in (I18ll9p . This leads to the 
following definition of the resonance curves: 



a^- 



{na + mj3y 



= /?^- 



2 [na + m/3)'^ 



= 0, 



0.2 




Fig. 5 Resonance curves in the a, (3 plane for the system H24 
with |n|, \m\ < 2. 



a 



±mj3 
2±n' 



P 



2±TO 



(23) 



for n,m = —00, . . . , cx), unless n{m) = ±2 in which case 
P{a) = 0. These are precisely the straight lines in the 
a, P plane where the frequencies are in ratio. 

As an illustration of how useful this technique can be 
we can examine the following problem, closely related 

to (uni. 



X +[a + e cos{pt)\x + e cos{at)y — 
y + e cos{pt)x + [/3 + e cos{at)]y = 0. 



(24) 



The time dependent terms here are of the same form as 
in (|16p. and thus we write solutions on the transition 
curves as in p7p . Now the recurrence relations, in the 
limit of e -^ 0, yield the following resonance curves: 



{na + mpy 



= [P 



[na + mpy 



= 0, (25) 



(22) 



for n,7Ti ~ —00,. ..,00 (we see that in this case there 
is no issue with convergence of the infinite Hill's de- 
terminants, as the driving frequencies do not equal the 
natural frequencies). The resonance curves for the first 
few values of n, to are given in Figure [5l which can be 
compared with numerically generated stability charts 
(using methods already discussed in Section 2) for non- 
zero e as shown in Figure [6] We see that the resonance 
curves as predicted by (|25p give a good approximation 
to the location of the bands of instability in parameter 
space, with very little effort. 
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0.25 
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Fig. 6 Stability diagram for system H24|l with e = 0.1 and the a, /3 unit square divided into an 800 x 800 grid. White represents 
points where the norm of the maximum multiplier is greater than 1.025, and black where the norm is less than 1.025. We see 
the strong correspondence between bands of instability in this plot and the resonance curves in Figure [5] 



In fact, we may generalize this technique to arrive 
at the fohowing proposition: 



Proof: On transition curves, the solution has the same 
quasiperiodic structure as Pi (t) and so we can write it 
as a Fourier series of the form 



Proposition: For the system 
^ + iPo + ePi{t))X = 



(26) 



X{t)=j2^ncxp 



n 



it 



(n.uj) 



with X e K", Po ^ diag{Xl, Xl, . . . , Xl) and Pi(t) of 
the form 



where An is an n-dimensional vector of Fourier coeffi- 
cients. Subbing this into (l^5|) we find 



Pi{t) 



'cos(ajii) . . . cos(w„i) ~ 



^cos(ajit) . . . cos(a;„i) , 



(n.o;)^ 



E^ 



nexp 



n 



it ■ 

— (n.u)) 



bands of instability will grow with e from resonance 
curves given by 



{Po + sPi{t))J2Anexp 



n 



it 



{n.u}) 



xi 



{n.Ljy 



0, j = l,...,n, 



where n is a vector of length n whose elements are all 
integers in the range — oo, . . . , oo, that is n S Z". 



The term in e in this equation is complicated and 
would involve shifting indices were we looking for the 
infinite Hill's determinant to determine the transition 
curves; however we seek only the resonance curves which 
are found by letting £ = 0. Canceling the exponentials 
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leaves 



E 

n 



{n.Lj) 



PnA 



o^n 



and since Pq = diag{Xf) for i 
with the n equations 

(n.o;)^ 



, n this leaves us 



A 



0, neZ" 



Now we may simply plot these curves for some elements 
of n as in Figure [5] 

We note this method predicts the resonance points 
for the Mathieu equation ([T]) to be (5— n^/4 = 0, and the 
resonance curves for the quasiperiodic Mathieu equa- 
tion © to be (5 - (n + muj)'^/A = (as in fill), with 
n,m & Z. 



4 The multiple scales method near the 1:1 
resonance 

Based on the previous section we see that the method 
of harmonic balance is unable to describe behavior near 
the 1 : 1 resonance, that is close to a = p. We use 
instead the method of multiple scales, which is a per- 
turbative technique that allows us to derive analytic 
representations near resonant values of a and /3. 

The typical way to proceed is to expand a, /?, x and 
y in powers of e, then define different time scales and 
separate out the behaviour on slower time scales. We 
could let 



ao + eai 



e^a2, 



/3 



ao + £/3i + 6^132 



so a and /3 only differ at first order. Then, the dif- 
ferential equations would contain terms involving a^t -\- 
aiet + a2e'^t which suggests we define three time scales, 

i, T^et, T = e^t 

(T here is not to be confused with the fundamental 
period T from Section 2). Note the e^ time scale is 
necessary as there are no secular terms at first order. 
Now we may expand x and y in powers of s where each 
term is a function of these three times, i.e. 



xo{t, r, T) -I- exi {t, r, T) -I- e^X2 {t, t,T) + . 



(27) 



and similarly for y. Subbing into the differential equa- 
tions and collecting powers of e we find a set of equa- 
tions with inhomogeneous terms involving partial deriva- 
tives of lower order solutions. Removing secular terms 



gives differential equations in the slow-time dependence 
of low order terms. However, we find that the crucial 
system of equations only describes flow on the 'very' 
slow time T, and is independent of ai and /3i (see the 
Appendix for details). 

Therefore, it is sufficient to define 



a 



a. 



I3^a + e'^l32. 



(28) 



Now there is only one slow time, T = e^t, however we 
are still expanding x and y in powers of e (rather than 
e^). To second order we find the following system of 
equations: 



e° -.-g^+a^xo^O, 

9^yo , 2 






+ a yo = 0, 



e" :-;tt7^ — I- a^xi = —xq cos{at + (i2T) — yn cos(ai), 
+ a^yi = ~yo cos{at + /32T) ~ yo coa{at), 



d'yi 



9t2 
2 d'^X2 2 o^^3:^o 



d^y2 



a X2 — —2 



+ a 2/2 



dtdT 
yiCos{at), 

d^yo 



- xi cos{at + /32T) 
yi cos{at + P2T) 



dtdT 

-yiCOs{at) -2a^22/o, 



where xq — xo{t,T) and so on. Solving one at a time, 
we find at zero order: 



xo = A{T) cos(ai) + B{T) sin(ai), 
yo = C{T) cos(ai) + D{T) am{at). 



(29a) 
(29b) 



There are no secular terms on the right hand side at 
first order so we find 



xi ^ a{T) cos{at)+b{T) sin(ai)- 



_A 
6a 



2 cos(2ai+/32T)- 



and a similar expression for j/i. Now when we look at 
second order, we see there are secular terms in the right 
hand side (that is, terms involving cos{at) and sin(ai)). 
Setting these terms to zero results in four linear dif- 
ferential equations in the slow time dependence of xq 
and yo, that is the four functions A{T), B{T),C{T) and 
D{T) defined in ^. Letting Y = {A, B, C, D) we may 
write these four equations as dY/dT = Q{T)Y where 



/ -sin(;32T) + 3sin(2/32r) -2 + cos(/32r) + 3cos(2/32r) 
2-1-5 cos(/32T) + 3 cos(2/32r) -lsin{fi2T) - 3 sin(2;32T) 
- sin(;32T) + 3 sin(2/32T) -2 -I- cos(/32r) + 3 cos(2/32r) 
V2 + 5cos(;32T) + 3cos(2/32r) -7sin(/32r) - 3sin(2;32T) 5- 



7sin(/32T) 
5 + 5 cos(/32T) 

7sin(/32T) 
24^3/32 + 5 cos(/32T) 



1 + cos(/32T) \ 

sHhT) 
24a3/32 + cos(/32r) 
sin(/32r) / 



Stability of a 2-diraensional Mathieu-type system with quasiperiodic coefficients 



11 



10 




0.05 



0.1 



0.15 



0.2 



Fig. 7 Tlie largest multiplier of system H30|l for /i small. 

The stability of this system will determine the sta- 
bility of the original system. 

This system is still Louiville {Tr{Q) = 0) and t- 
invariant (as defined in (|TT|) with G = diag{l, —1, 1, — 1)). 
A nice simplification happens if we define a new time 
coordinate t* = P2T, and we may write the system in 
the form 

— = [A + ^(Ao + Ai cosit*) + A2 cos(2t*) 
+Bi sm{t*) + B2 sin(2t*))] X 



(30) 



where fi — l/(24a'^/32) and the constant coefficient ma- 
trices are given by 



A = 



/ \ 



1 

V -1 0/ 



Ao = 



/ -2 1 \ 

2 5 

0-201 

V 2 5 0/ 



with Ai, A2, Bi, B2 in the Appendix. This is a one pa- 
rameter 27r-periodic system, with the n parameter mul- 
tiplying terms very much like a Fourier series. It is 
straightforward to analyze the stability of this system 
using the methods already outlined, and we show the 
maximum multiplier in Figure [T] 

The appearance of a small region of stability can 
be understood in the following way. In system (I5(I)) . 
when /i = the system has natural frequency 1. As 
we switch on /x, the driving frequency is also 1 and res- 
onance causes the system to be unstable. However, as fx 
grows, the natural frequencies are now the eigenvalues 
oi A + iiAq , which become sufficiently distinct from the 
driving frequency for stability to resume. This is tran- 
sient however and as the 'perturbation' /u grows the sys- 
tem becomes unstable again. This window of stability 
occurs for 

{fi- ^) 0.18106 <fi< 0.189 (= ^+), 

which tells us the values of P2 where stability occurs 
and this in turn gives us the following band of stability 



CQ. 0.5 - 



0.25 




0.25 



0.5 



0.75 



Fig. 8 The multiple scales prediction for a band of stabil- 
ity overlaid onto the stability chart generated numerically in 
Section 2. 



in (a, /3) parameter space 



a± e' 



( 1 



V24a3/i4 



< /3 < a ± e' 



/ 1 



V24a3//_ 



This result is in excellent agreement with the numerical 
results generated earlier, as shown in Figure |S1 at least 
for values of a, /3 over 0.5. As discussed previously, this 
is because small values of a,/3 are equivalent to large 
values of e and the multiple scales method would not 
be expected to work well in this regime. Note that the 
symmetric nature of the system is recovered by the ex- 
pansion 



/3- 



0.2-, 



(31) 



analogous to (|28p. which results in the same stability 
region as in Figure [5] reflected about the line a = /?. 



5 Conclusions 

In this paper we have begun to analyze the rich dynami- 
cal structure of two dimensional systems with quasiperi- 
odic coefficients. Through a combination of numerical 
and perturbative techniques we are able to predict and 
determine the transition curves in parameter space which 
demarcate the regions leading to stable and unstable 
solutions. In particular, we find that using the method 
of harmonic balance we can quickly approximate the 
resonance curves in parameter space from which bands 
of instability emanate, with the caveat that the band 
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of instability may spread out on either side of the res- 
onance curve or move to its side. The muhiple scales 
method has proven to be particularly informative, if a 
little cumbersome, for examining the vicinity of specific 
resonances. 

There are three obvious issues which need to be un- 
derstood further with regards the particular system un- 
der consideration in this work. The first is the 'arcs' 
of instability which are commented on in §2 and can 
be seen in Figure [51 These arcs do not arise in the 
harmonic balance analysis, nor do they correspond to 
a linear relationship between the frequencies and thus 
the method of multiple scales would be inappropriate. 
The exact mechanism which leads to the development 
of these arcs is not understood, and is worthy of further 
investigation. 

The second issue is the possible resonances between 
first order terms in the expansion of a, /3 as described 
in the Appendix. Perhaps these could explain why the 
multiple scales method predicts a continuous band of 
stability near the 1:1 resonance, whereas numerically 
we observe a broken series of islands of stability (see 
Figure [2] and [8]) . Finally, all bands of instability visible 
in Figure [2] taper as we increase a and /3, as predicted 
by the harmonic balance method, with the exception 
of the large region of instability surrounding the 1:1 
resonance. While the multiple scales analysis could be 
interpreted as telling us that the true band of instability 
around a = /3 is in fact very narrow and exists between 
the predicted bands of stability, as shown in Figure |S1 
it still leaves an open question as to the origin of this 
large region of instability. 

There are numerous extensions to the system in 
question, however it would be wise to keep the number 
of free parameters small to enable an informative anal- 
ysis. For example, a general form could consist of two 
natural frequencies and four driving frequencies. While 
this would undoubtedly lead to exceedingly rich sta- 
bility properties, six free parameters would be difficult 
to present in a clear fashion. A 'detuning' parameter 
in the driving frequencies would instead lead to three 
parameters and would provide an interesting extension. 
There are of course other extensions involving nonlinear 
terms, more dimensions and so on. Finally, there are 
numerous applications involving quasiperiodic oscilla- 
tions about fixed points whose analysis would benefit 
from the present work. 



6 Appendix 

As is typical in multiple scales analysis, the equations 
at high orders are quite lengthy and involve very many 



terms; as such we will only give an overview to justify 
ignoring first order terms in a, /3 to arrive at (|28p . 
Let us write 

a — ao + eai + e^Q;2, /3 = ao + eA + ^'^^2 

and expand x,y in the three time scales t, t and T, as in 
(P7)) . The second time derivative becomes (up to second 
order in e) 



d^ 



d^ \ 



d^ d^ 



dt^ dt^ '^y^dtdTj^^ X^'dtdT ' arV 

leading to the following system of equations in x: 



^ a;o 2 n 

1 d'^xi ^d^xo 2 

s ■-Q^ + ^Q^+'^o^i + ^'^oaiXo 

+ xo cos(aot + /?iT + I32T) 

+ yo cos(Q;ot + aiT + a2T) = 

9 d^X2 ^d^Xi ^d^Xn d'^Xn 

P^ ■ 1 -|_ 2 - + 2 - -I - 

■ dt^ dtdr dtdT dt^ 

+ aoa;2 -I- 2aoaiXi + {al + 2aoa2)xo 

+ xi cos{aot + /?iT + I32T) 

+ Hi cos(Q;ot + aiT + a2T) = 

For the corresponding system in y, simply swap x 
with y, and ai,a2 with /3i,/?2. We examine the equa- 
tions one order at a time. At zero order, the solutions 
are 

xq — a{T, T) cos(Q;ot) -f 6(t, T) sin(aot), 
yo — c{t, T) cos(aoi) -f d[T, T) sin(aoi). 

At first order, which we write as 

d'^xi 



a^xi = inhomogeneous terms. 



the right hand side includes 

2ao sin(Q:oi) I ^ aih j + 2ao cos(aot) I —-z aia 



in the x-equation, and similarly for c, d in the y-equation. 
Setting these secular terms equal to zero gives us the 
r-dependence at zero order, i.e. 

a{T, T) = A{T) cos(aiT) + A' {T) sin(aiT), 
5(r, T) = B{T) cos(air) + B' {T) sin(aiT), 

and similarly for c, d. The solution at first order will 
then involve the homogeneous solution 

xi = P{t, T) cos(aot) + Q{t, T) sin(aot), (32a) 

2/1 = R{t, T) cos{aot) + S{t, T) sin(aoi), (32b) 
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and the inhomogeneous solution, made up of sin and 
cos terms with arguments 

(ao ± ao)t + (ai//3i)T + {a2lP2)T 

which we wih not hst in full (by ai/ Pi we mean linear 
combinations of the form fciQ^ + ^2/3^ with /c S Z and 

|fcl| + |fc2|<3). 

At second order, which again we write as 



9^ 



OioX2 — inhomogeneous terms. 



we use the zero and first order solutions to calculate 
the right hand side and look for secular terms, that 
is terms multiplying cos(Q;ot), sin(Q;ot). The full expres- 
sion for the right hand side is very lengthy and compli- 
cated, and involves multiplying many cos and sin terms 
of different arguments together, however we are only 
interested in those with ag in the i-dependence and so 
can ignore those containing 2ao or Sag. This gives us 
four equations: the coefhcients of cos(aoi),sin(aot) in 
the x-equation, and the same for the y-equation, which 
must be set equal to zero to prevent resonance in X2,y2- 
These four equations are of the form 



do 
'2ao^ - 2aoaiP + h{t, r, T) - 0, 

OT 

dP 

2aoT^ 2aoaiQ + f2{t,T,T) = 0, 

OT 



(33a) 
(33b) 



and again in R,S and fs, fi with ai replaced by /?i. 
These define the r-dependence oi P, Q , R, S (from the 
first order solutions xi, yi in ([35])), and we see these will 
be oscillatory if we remove the cos(air),sin(air) terms 
in /i, /2, and the cos(/3iT),sin(/3iT) terms in /a, f^. This 
leads to eight equations involving the eight functions 
describing the T-dependence of the zero order solutions, 
that is 

A,A',B,B',C,C\D,D' 

and their derivatives with respect to T. We may write 
this as a first order system dZ/dT = A{T)Z, where A 
has block form 



A{T) 



Ai A2{T) 
A^{T) Ai 



with A being Louiville (vanishing trace) and i- invariant, 
as defined in (ITTt with 

G = diag{l, -1, -1, 1, 1, -1, -1, 1). 

The stability of our original system now hangs on the 
stability of this one, and crucially the parameters ai, (3i 
do not appear in the matrix A. Thus the difference be- 
tween a and /3 at first order is irrelevant and we may 
set ai = f3i — 0. 



We note that there are a number of possible reso- 
nances in the first order terms in a, /3 which we have 
ignored. The fi functions described above contain terms 
like, for example, 



cos(2aiT + /3iT + a2T - ^T). 



(34) 



In general, this will not be secular in equations p3p . If 
however it was the case that /3i = —ai then this would 
be secular. There are six combinations which lead to 
additional secular terms, they are 



ai 



±/3i, Qi ~ ±3/3i, 3ai = ±/3i. 



Each would need to be considered in turn, and would 
lead to a different coefficient matrix A given above. We 
will not pursue this issue further in the present work. 
Finally, we give the coefficient matrices from (PH)) . 



Ai 



5 



V 5 

/ 



Bi 



V 



1 1 \ 

5 

1 1 
5 0/ 

7 0\ 
-7 1 

7 
-7 1/ 







-1 



A. 



,B. 



/ 3 \ 

3 

3 

V 3 / 

/3 0\ 
0-300 
3 

VO -3 0/ 
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